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UNAMIT - A ONE-DIMENSIONAL 4tt SPHERICAL MULTILAYER 
REACTOR-SHIELD-WEIGHT OPTIMIZATION CODE 

by Eugene S. Troubetzkoy, * and Millard L. Wohl 
Lewis Research Center 

SUMMARY 

The computer code UNAMIT allows the computation of total neutron and gamma- ray 
dose rate at the surface of a multilayer spherical shield surrounding a spherical reactor 
model whose leakage is known. It computes this dose rate subject to a minimum total 
shield weight constraint. The code utilizes an exponential attenuation - secondary- 
gamma- ray production model with a transport-determined dose rate and artificial source 
term to generate the variation of dose rate through the shield and obtain its value at the 
shield surface. 

The user must specify a set of parameters characterizing the geometry and com- 
position of the reactor- shield configuration. These are (1) the reactor outer radius, 

(2) the dose rate convergence criterion, (3) minimum values for all shield layer radii, 

(4) maximum values for all shield layer radii, (5) step sizes for shield layer radii, and 
(6) shield layer material densities. 

A set of parameters characterizing the materials specified in the shield layers is 
also required. These must be specified by the user and are for each shield layer: 

(1) attenuation parameters for primary gamma rays, p, (2) attenuation parameters for 
neutrons, y, (3) attenuation parameters for secondary gamma rays, X, and (4) secondary 
gamma-ray production parameters, z. These parameters depend not only on the shield- 
ing materials in question, but also on an assumed preliminary shield configuration. 

The model assumes exponential attenuation and secondary production in each shield 
layer, and is designed to match the neutron dose rate at the shield surface as determined 
by a transport calculation. The parameters used by the calculational model are obtained 
from preliminary transport analyses so that the calculational model gives a reasonable 
approximation of the attenuation and secondary production of radiation. 

The UNAMIT code adjusts the thicknesses of the specified shield layers by finite in- 
crements within layer boundary limits prescribed by the user. It generates, by a step- 
ping search procedure, the minimum weight shield configuration satisfying the shield sur- 
face dose rate constraint specified by the user. The minimum weight configuration obtain- 
ed lies within the set of possible configurations defined by the specified shield layer mate- 
rials and their order, and the minimum and maximum prescribed shield layer radii. 

* Formerly at United Nuclear Corporation, Elmsford, New York; presently at 
Columbia University, New York, New York. 



INTRODUCTION 


The design of in unit nuclear reactor shields to permit the achievement of low 
dose rates in the near vicinity of the powerplant is an involved, time-consuming proc- 
ess. Detailed neutron and gamma- ray transport calculations are an important part of 
the overall design procedure. If the constraint of minimum unit shield weight is added 
to the required dose rate constraint, the shield design becomes more difficult. 

Weight- optimization of even first-order unit shield configurations is a highly in- 
volved procedure if iterative neutron and gamma-ray transport analyses are done. 

These are very time consuming and the probability of error is high. 

An alternative to iterative transport procedures is to use a highly simplified com- 
putational model. Such models were first systematically employed in the Military Com- 
pact Reactor (MCR) Program where system constraints necessitated minimum shield 
weight. During the course of this program, a series of simplified shield weight- 
optimization codes was developed. 

These codes used a model starting with a reactor leakage source term and an ex- 
ponential attenuation - secondary-production description. Thicknesses of prescribed 
ordered shield layers were systematically varied to achieve minimum shield weight. 

The latest variant of these codes is called UNAMIT. This code was used exten- 
sively in the nuclear airplane studies being conducted by the NASA Lewis Research Cen- 
ter. The code accepts an ordered set of shield materials, nominal inner and outer layer 
radii, and radial step size increments. It then generates the minimum-weight layered 
unit shield configuration consistent with the prescribed shield surface dose rate. 

In order to generate a lowest-weight shield configuration, UNAMIT requires atten- 
uation and production parameters characteristic of the shield materials and their order 
in the shield. The required parameters are attenuation parameters for primary radia- 
tion; attenuation parameters for secondary-producing radiation; attenuation parameters 
for secondary radiation; and secondary-production parameters. These parameters de- 
pend upon an initial assumed shield configuration and are generated from preliminary 
transport analyses for such a configuration. 

This report will discuss the operation of the UNAMIT code and the pertinent theory 
involved. Listings of the code in both Fortran IV and Fortran 63 are presented. In ad- 
dition, full input instructions and a sample problem output listing are included. 


ANALYSIS 

The code UNAMIT computes the total neutron and gamma-ray dose rate at the sur- 
face of a multilayer-unit shield as described below in equation (1). It weight-optimizes 
the shield subject to a fixed dose rate constraint at the shield surface. The following 
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section describes the computational method, including the determination of primary and 
secondary attenuation and production parameters. 


Calculational Procedure 


The calculational model uses an arbitrary number of neutron groups, each of which 
produces secondary gamma rays. An arbitrary number of primary gamma-ray groups 
is also specified. The dose rate at the shield surface is the sum of the neutron dose 
rate and the primary and secondary gamma ray dose-rate contributions. 

For a spherical shield consisting of N layers of distinct materials between radii 
r 0 , r 1? r 2 , . . r n , the dose rate is given by 


D = 
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where 

i outer summation index denoting region with inner and outer radii r^_^ and r i 

i inner summation index for regions internal to r n region, which affect produc- 

J th 

tion of secondary gamma rays in i m region 

K energy group summation index for secondary-producing neutrons and primary 
gamma rays 

GS^ source intensity for neutrons emitted at Tq which are responsible for secondary 
gamma rays produced in region i; for each region denoted by i, K runs from 
1 to NGS.. 

attenuation parameter for K th neutron group in region j for fixed i; it is de- 
fined only for j < i 

Xtr attenuation parameter for gamma rays produced in region i by K neutron 
group; it is defined only for j > i 

Z? secondary gamma-ray production parameter, photons/(source neutron-cm) for 

1 ru 

K in neutron group in region i 
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ir 

GP source intensity for primary neutrons and gamma rays emitted at r^ which con- 
tribute to external dose 

NGS^ number of secondary gamma- ray producing groups in region i 
TC th 

PP attenuation parameter of K primary gamma- ray group in region i 
NGP number of primary radiation groups 
N number of layers of distinct materials 

The shield weight is then given by 



(2) 


where is the i^ 1 region material density. 

Discussion of equation (1) . - The first sequence of terms in equation (1), within the 
double summation, deal with the production and subsequent attenuation of secondary 
gamma radiation. 

GSx 1 is the source intensity for reactor leakage neutrons which are responsible for 
secondary gamma rays produced in region i. The exponential term 


exp 


i-1 

-I 


Y ji (r j - Ti> 


describes the attenuation of these neutrons in the outward radial sense through shield 
regions up to, but not including region i, where secondary gamma rays are produced. 


The exponential term exp 


Y^(r 


i-1 


describes the leakage neutron attenua- 


K 


tion through region i, in which secondary gamma rays are produced. The term 

is the production strength for secondary gamma rays produced in region i. 

r k" I 

The exponential term exp - Xf:(r. - r) describes the attenuation of secondary 
gamma rays produced in region i through the remainder of region i; the product of 
doubly subscripted terms (subscripts ii) must be included in the integral over region 
because the secondary production may take place anywhere in this region. 

The exponential term 


exp 


~ N 
- £ 


j=i+l 


- r M> 
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describes the attenuation of secondary gamma rays produced in region i through shield 
regions external to i. 

The last term in equation (1) 


GP K exp I 


N 


-I 

i=l 


P^( r i - r 


i-l' 


describes the source strength and attenuation of primary radiation, both neutrons and 
gamma rays. 

Discussion of attenuation and secondary-production parameters . - The model de- 
scribed by equation (1) assumes that a given shielding material is described by a num- 
ber of parameters. These are the attenuation parameters for primary radiation P, the 
attenuation parameters for secondary producing radiation Y, attenuation parameters 
for the secondary radiation X, and secondary production parameters Z. These param- 
eters depend not only on the shielding materials in question but also on the shield con- 
figuration. It is, therefore, important to obtain the model parameters from transport 
calculations so that the model gives a reasonable approximation of the attenuation and 
secondary production of radiation. 

The methods of determining the primary and secondary radiation attenuation and 
production parameters will be discussed using a two-layer shield model consisting of 
an inner layer of high-atomic number gamma-ray shielding material followed by a layer 
of hydrogenous neutron shield material. 

Primary radiation . - 

Neutron dose rate contribution: In order to determine the neutron attenuation pa- 
rameters, let us consider a single group of neutrons responsible for the neutron dose 
rate at the shield surface. This will be the first category of primary radiation. A 
transport (Monte Carlo or Sn) calculation of the total primary neutron dose rate as a 
function of r, the distance from the reactor center, might yield the following plot. 
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The dotted lines are straight lines having the same slopes as the solid curves at r, and 

11 1 
r£. The slopes are P^ and Pg respectively, and serve as the attenuation parameters 

for primary neutrons in the model. The ordinate intercept, In Gp\ is obtained by ap- 
pending to the dotted line segment de a segment cd parallel to ab. 

Given all the attenuation parameters, P^, the source term GP* is uniquely deter- 
mined as indicated and is described, for the configuration of sketch (a), in the equation 


GP 1 e 


- p i (r r r o)J L- p 2 (r 2 ' r i ) 


= D 


1 2 

where D = 47rr x the primary neutron dose rate at the shield surface. 

Primary gamma-ray dose rate contribution: The primary gamma-ray contribution 

o 

is designed by the index K = 2 in equation (1). The parameters P^ for primary 
gamma-ray attenuation can be determined in a manner similar to the method just dis- 
cussed for obtaining primary neutron attenuation parameters. An alternative would be 

to use the mass absorption coefficient, for an appropriate gamma-ray energy, for the 

o 

material in region i as P. . 

Secondary radiation . - In order to determine the secondary gamma-ray attenuation 
parameters, let us assume that a single group of neutrons is responsible for the produc- 
tion of a single group of secondary gamma rays. 

All of the secondary gamma-ray attenuation coefficients At: can be assumed to be 
given by the mass absorption coefficients in the material for a representative energy of 
gamma rays produced in region i. 

From the results of a transport calculation, the representative energy of neutrons 
responsible for secondary production in region i can be determined. By fitting the 
model of the attenuation of the neutron group in a manner similar to the one described 
under "Primary radiation, " the secondary- producing neutron attenuation parameters 
Yf: in region j may be determined. 

J 1 pr 

A special procedure can be applied to determine the attenuation coefficient of 

neutrons in region i (either region shown in sketch (b)) responsible for secondary pro- 
duction in the same region i. A transport calculation can give D^, the contribution to 
the shield surface dose rate due to secondary production at all positions r in the shield. 
By plotting ln(47rr x this contribution) against r, a curve which should be fit by a 
piecewise linear curve (according to our model) is obtained as follows: 
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y-Shield region Neutron shield region 

tb) 


In each shield region, the slope of the straight line fit is obtained from 
Y-| = slope + X^.. Having obtained X- from the mass absorption coefficients in the 
manner already prescribed, one can obtain Y. . from the slope. 

Finally, the combination of parameters GS? • Z ? can be determined by back- 
extrapolating the total secondary dose rate due to region i as shown in sketch (c). 



In sketch (c), . is the secondary dose rate from region i. As in sketch (a), cd is 

parallel to ab. 

As has been mentioned throughout this discussion, useful application of the 
attenuation-production model is dependent on an accurate transport calculation for a 
shield similar to the set of shields to be studied with the UNAMIT code. It is recom- 
mended that the accurate calculation be performed by either Monte Carlo or Sn methods 
employing carefully prepared cross sections and/or transfer matrices. 
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ASSUMPTIONS AND LIMITATIONS OF UNAMIT 


As has been mentioned, UNAMIT utilizes an exponential attenuation - secondary- 
production model. The input parameters characterizing the initial assumed shield layer 
configuration are derived from transport analyses for an assumed physical shield layer 
configuration. If, using this assumed initial configuration, UNAMIT cannot obtain the 
surface dose rate constraint prescribed by the user, it may be necessary to assume a 
second initial configuration from which new attenuation and production parameters must 
be obtained. Shield calculations made to date have never required a third initial con- 
figuration. 

As an adjunct to the previous statements, the user should be careful not to assume 
that attenuation and production parameters which were successful in generating optimum 
shield weight configurations in a given reactor power regime would be successful in an- 
other regime. Changing reactor power regimes is equivalent to changing attenuation 
regimes, and generally new parameters are required when this is done. 


PREPARATION OF INPUT DATA 

The user specifies minimum and maximum values of all radii and step sizes for 
each radius. UNAMIT systematically varies each radius within the prescribed range, 
adjusting them to give the prescribed dose. The shield weights for all configurations 
within the variation considered are compared with one another, and the configuration 
yielding the lightest shield is retained. 

The input description is as follows: 

INPUT 

PART A 1st Card: Shield configuration and specifications 

Title (any 80 Hollerith characters) 

2nd Card: N, number of regions 

Tq, core radius 

D, dose specification - left blank for UNAMIT FORMAT (I 12, 
2E12. 5) 

The next_group of cards must be repeated for regions i = 1, 2, 3, . . . , N: 

1st Card: i, region number 

NGS^, number of distinct secondary production processes 
occurring in region i 

p., density of material in i^ 1 region (g/cm^) FORMAT 
(215, E 12. 5) 
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2nd Card: 

3rd Card: 

3 + j th Card: 

After the N groups 
primary radiation: 

1st Card: 

2nd Card: 

3rd Card: 

3 + j th Card: 
PART B 1st Card: 


GS?, source term for group K in region i 

Zf , secondary production parameter for group K in region i 
K = 1, 2, . . . , NGS-, five pairs of numbers to a card 
FORMAT (10E8. 4) 

Y^h, neutron attenuation parameter 

K* 

x£ ., gamma-ray attenuation parameter (blank if i > l) 

K = 1, 2, . . . , NGS. , five pairs per card. Note: for each 
K value, a set of values of Y and X is specified for 
i + 1 < l < N. 

FORMAT (10E8. 4) 

YjT, neutron attenuation param«X~- (blank if j > i) 

Xjj, gamma-ray attenuation parameter (blank if j < i) 

K = 1, 2, . . . , NGSp five pairs per card) 

FORMAT (10E8. 4) 

of cards have been entered, enter the following information on 


NGP, number of primary groups 
FORMAT (112) 

GP*^, primary source terms, K = 1, NGP (10E8.4) 

TC th 

P[% attenuation parameter of K primary radiation in 

region i 

K = 1, 2, 3, . . . , NGP 
FORMAT (10E8. 4) 


Pt , attenuation parameter of K primary radiation in 
region j. K = 1, 2, 3, . . . , NGP 
FORMAT (10E8.4) 
j runs from 1 to i 


EPS, convergence criterion 

NITMAX, maximum number of iterations permitted 
ISW; 0, no intermediate answers 

1, answers for all shields calculated 
FORMAT (E12. 4, 216) 
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2nd Card: RMI(l) RMT(2), . . . , minimum values for all radii 

FORMAT (6E12. 4) 

3rd Card: RMS(l), RMA(2), . . maximum values for all radii 

FORMAT (6E12. 4) 

This can be followed by PART A and PART B for as many shields as desired. 


OUTPUT 


The input information is printed out. 

If ISW = 1, radii and weights of all shields calculated are printed out. 
The lowest weight shield obtained is edited. 


SAMPLE PROBLEM 

Monte Carlo Analysis of Reactor-Shield Configuration 


A spherical homogeneous 375- megawatt reactor with a layered uranium-water sphe- 
rical shell 47 r unit shield was analyzed with the SANE-2M and SAGE 4 Monte Carlo codes. 
SANE-2M is used to perform neutron transport calculations and generate secondary gam- 
ma rays. SAGE 4 is used to perform primary and secondary gamma ray transport calcu- 
lations. 

The configuration analyzed with SANE-2M and SAGE 4 is indicated in table I. 


TABLE I. - CONFIGURATION FOR MONTE CARLO REACTOR- SHIELD ANALYSIS 


Physical region 

Type 

Outer radius, 
cm 

Physical region 

Type 

Outer radius, 
cm 

1 

Core 

82. 38 

14 

h 2 ° 

150 

2 

Ba reflector 

90 

15 

u 

152 

3 

U 

93 

16 

h 2 ° 

160 

4 

H 2° 

100 

17 

u 

162 

5 

u 

103 

18 

H 2 ° 

170 

6 

h 2 ° 

110 

19 

u 

172 

7 

u 

113 

20 

H 2 ° 

190 

8 

H 2 ° 

120 

21 

u 

192 

9 

u 

123 

22 

H 2 ° 

210 

10 

H 2 ° 

130 

23 

u 

212 

11 

u 

132 

24 

H 2° 

230 

12 

H 2 ° 

140 

25 

h 2 ° 

320 

13 

u 

142 

26 

Air 

1000 
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TABLE II. - REACTOR- SHIELD MATERIAL 


COMPOSITIONS 


Component 

Element 

Atom density 
atoms /cm ^ 

Core 

Hydrogen 
Oxygen 
Aluminum 
Zirconium 
Uranium- 2 35 

1. 976X10 22 
1. 184X10 22 
5. 12X10 21 
1.744X10 22 
9.79X10 20 

Reflector 

Beryllium 

1. 20X10 23 

Depleted uranium 

Uranium- 238 

4. 82X10 22 

Water (borated) 

Hydrogen 
Oxygen 
Bor on- 10 
Boron- 11 

6.45X10 22 

99 

3. 37xl0 z ^ 
1.73X10 20 
7. 85X10 20 

Air 

Nitrogen 

Oxygen 

4. 17X10 19 
9.0X10 18 


The material compositions for the configuration specified in table I are listed in table II. 

A Monte Carlo analysis of this reactor shield configuration yielded a dose rate of 
0. 0224 millirem per hour at a detector point in air 30 feet from the reactor center. 


Optimum Shield Weight Determination With UNAMIT Code 

The Monte Carlo analysis was performed with specific uranium- water configurations 
represented by discrete uranium layers separated by water. The uranium layers were 
of 2- to 3- centimeter thickness to allow the benefit of self- shielding effects. The thick- 
ness of the water layers was chosen to yield effective mixture densities in the range of 
3 to 6. 5 grams per cubic centimeter. 

In performing shield optimization analysis with UNAMIT, three homogenized 
uranium-water (borated) mixtures were considered in the shield in addition to an outer 
layer of borated water. The mixtures were 

(1) 30 Percent depleted uranium plus 70 percent borated water, giving a density of 

6-. 4 grams per cubic centimeter 

(2) 20 Percent depleted uranium plus 80 percent borated water, giving a density of 

4. 6 grams per cubic centimeter 

(3) 10 Percent depleted uranium plus 90 percent borated water, giving a density of 

2. 8 grams per cubic centimeter 
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Mixtures 1, 2, and 3 correspond to homogenized combinations of physical regions 3 to 
10, 11 to 18, and 19 to 24, respectively, in table I. 

The weight- optimized configuration generated by UNAMIT, for a dose rate con- 
straint of 0.025 millirem per hour 30 feet (9. 14 m) from the reactor center, is shown 
schematically in sketch (d). 


Core and jy Mix 


reflector / / ture 1 


Mix- 

Iture M ,xture 3 


Radius, cm 90 93 


2 

130 155 


Water (borated) 


215 


322 


id) 


The weight of this optimized configuration is 543 000 pounds (247 000 kg) for the pre- 
scribed 30-foot (9. 14- m) dose rate of 0. 025 millirem per hour. This is very close to 
the Monte Carlo-determined dose rate of 0. 0224 millirem per hour at the same detector 
location. Details of the computational results are listed in the following output listing. 
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FORTRAN IV COMPUTER CODE LISTING 


C PROGRAM UNAMiT 

COMMON /MA X/N , Z ( 2) »RC19) ,X<2) ,CON,W 

DIMENSION RMI (20) ,RMA (20) ,DER(20) , ROPT (20) , RNG( 20 ) 

100 CALL IN IT AL 

R EAD( 5 9 10 1 ) EPSfNITMAXtlSW 
READ( 5,102) (RMI ( I ) ,1 = 1,N) 

WRI TEU, 1C2) (RMI < I ) ,1 =1,N) 

REAO ( 5, i 02) ( RMA( I ) ,1 =1 ,N) 

WR I TE ( 6 , 1 02) (RMA( I ) ,1 = 1,N) 

READ ( 5, 102) (DER( I ) ,1 = 1,N) 

WR I TE ( Ey 1 02) ( DER ( I ) ,1 = I,N) 

101 FORMAT (E 12*4,216 ) 

102 FORMAT (6E12*4) 

WMAX*- UOE + 25 
NM 1=N- 1 

DO 1 1=1, NM1 
R ( I ) =R M I ( I ) 

1 RNG ( I ) =RM I ( N ) 

R ( 1 ) = R ( 1 )“DER ( 1 ) 

2 1 = 1 

3 R ( I )=R ( I ) +DER ( l) 

1 F ( R( I HRMA(I ) ) 4,4,15 
A I F( R( I )-R(I + l))5t 5»15 

5 RN1 =RNG( I ) 

R ( N )=RN 1 
NIT=NI7MAX 
CALL DR 

CON 1=CCN~ 1*0 

IF (ABS (CON1)-EP$)10*10,6 

6 RN 2= RNUDER(N) 

R ( N )=RN2 

CALL DR 
CON 2=CCN- 1 e 0 

IF ( ABS CCGN21-EPS) 10,10*7 
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? RN = RN1-C0N1* CRNI-RN2)/ICCN1h:ON2I 

RCN l=RN 
CALL OR 
CON =CG N- 1 « 0 

IFCABS (CON hEPSU0a0f8 

8 NIT=Ni1-l 
IriHll )99 9 99 t 9 

99 STOP 

9 RN 2=RN 1 
CON 2 = CCN i 
RN 1=RN 
CON 1=CCN 
GO TO 7 

1C RNGU ) =R ( N ) 

IF(RtN)-RCN-l) ) 2, 11 * 1 1 

11 CALL WR 

I F( ISW ) 12 » 12 v 1 1 1 
J1 1 WR I TE ( 6,1 12) W, (R( J) , J=1,N) 

112 FORMAT ( 10E12*4) 

12 I F ( W- In MAX ) 2, 2 , 13 

13 WMA X= W 

DO 14 J = 1,N 

14 ROP T( J )=R ( J ) 

GO TO 2 

15 R ( I )=RM I ( I ) 

IF ( I- 1 199 » 18, 16 

16 I F( R( I )~R ( I- 1 ) >17,18,18 

17 R( I )=R ( I- 1) 

18 1=1+1 

I Ft I-N >3, 19, 19 
19 DO 20 1=1,20 
2 C R ( I )=RQPTU) 

CALL EDIT 
GO TO 100 
END 


$ I B F 7 C IN I T UN DECK 


SUBROUTINE INI TAL 

COMMON /MA X/N, MNGP ,ALAM,R< 19 > , QZZZ , CSQ ,CGN , W 

DIMENSION R0< 19) »G S (1 0 ) ,Z ( 10) , PARI 10) ,61100) , PARP( 10 ) , P( ID 0, 19 ) , 
1G< 19), DOS! 100) 

DIMENSION TITLE(IC) 

999 FORMA T ( 1H 1 ) 

WRITEi £,9 99) 

ALPHA = (4.C*3« 14159} /( 453* 5924E+5*3. ) 

1000 FORMAT ( 10A6) 

READ (5, 1000) (TITLECI ) ,1=1,10) 

WR I TE ( 6 , 1 000 ) ( TI TLE (I) , I = 1 , 10 I 

1002 FORMAT ( 3HGN=I 12 , 7H RZER0=E12. 5 ,6H D8AR=E12e5) 

2002 FORMAT (I12,2E12«5) 

REAO ( 5 , 2C02) N,RZERC, DBAR 

WRITE( 6, 1002) N,RZER0,D8AR 

L =0 
M = 0 
KP = 0 


DO 17 1 = 1, N 

1003 FORMA T ( 4H0KR = I 5,4H NG=I5,4H RC=E12«5) 
READ ( 5 , 2 CC3 ) KR , NG , ROKR 

2C03 FOR MAT (21 5,E12»5) 

. ROC KR ) =ROKR 



KR,NG,RC!KR) 


WRITE ! 6 ? 1 C 03 ) 

I F(KR-i) 901,1,901 

901 WRITE ( C, 3901 ) 

3901 FORMA T ( 2 1 H REGION NUMBERS NOT CONSECUTIVE) 

STOP 

1 READ (5,1005) C G S ( K ) , l C K) , K= 1 , NG ) 

1004 FORMA T ( 5E 1 2 * 5 ) 

WRITE! 6,998) 

558 FORMAT! 4H GS) 

WR I TE ( 6,1004) (GS( K) ,K=1 , NG ) 

WRI TE ( 6,9 97) 

997 FORMA T ! 4H Z) 

WRITE! 6,1004) ( Z ( K ) ,K=1 ,NG) 

WRITE! 6,956) 

596 FORMA T ( 4H PAR ) 

N G= 2*N G 
DO 16 J = 1 , N 

READ ( 5,1005) ( PAR ( K ) , K=1 , NG ) 

ICO 5 FORMAT ( 10E 80 4 ) 

WRI TE ( 6, 1006) ( PAR ( K ) ,K=1,NG) 

ICO 6 FORMAT ( 10E12.5) 

IF(J-l) 902,2,4 

902 WR I TE ( 6,3902) 

3502 FORMAT ( 16H J LESS THAN ONE) 

STOP 

2 DO 3 K = 1, NG, 2 
K P=KP + 1 

ZZZ=EXP!RZERO*PAR!K) ) 

IF (ZZZ.GT.1.0E+27) GO TO 111 
IF! G S ( KP) oGTeloCE+10) GO TO 111 
2 GS ! KP )=GS(KP)*ZZZ 
GO TO 199 

111 GS! KP ) = i.OE+37 
199 CONTINUE 
KP = 0 

4 IF! J- I ) 5,7,12 

5 DO 6 K = 1, NG , 2 

6 PAR(K+1)=PAR(K) 

GO TO 5 

7 DO 8 K= 1, NG , 2 
K =K 

M =M + 1 
KP =KP + 1 

DEBUG K, A, (PAR! JD) , JD=1 , 10) 

A = -GS(KP)*Z(KP)/(PAR(K )-PAR(K +1)) 

DEBUG A 
G(M ) = A 
M =M + 1 

IF! I- 1 ) 90 2 , 5 03 , 5 04 

90 3 XXX=E XP t { R ZERO# ( PAR (K+ 1 ) -PAR ( K) ) ) ) 

IF .(XXX.GT.1.0E + 27) GO TO 777 
IF !A®GT.1«0E*1G) GO TO 777 
G! M )=-A*XXX 
GO TO 779 
777 G ( M )= 1 «CE+37 
779 CONTIN1E 
GO TO £ 

904 G! M )=-A 

8 CONTINUE 
K P - 0 

5 IF! J- 1) 9C2, 10,14 
1C DO 11 K= 1 ,NG 

11 PARP! K )=P AR! K i 
GO TO ] 6 

12 DO 13 K ~ 1 ? NG , 2 
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13 PAR(K)=PAR(K«-1) 

14 DO 15 K = 1 ? NG 
L =L + 1 

P(L?J— 1)=PARP(K)-PAR(K) 

15 P AR P ( K )=PAR*K > 

L=L-NG 

16 CONTINUE 

DO 1? K=1,NG 
L =L + 1 

17 P ( L ,N ) =PA R { K ) 

READ* 5,2002)NG 

ICO 7 FORMAT (4HCNG=I 1C) 

WR I TE ( 6, 10G7JNG 

READ( 5, 1CC5) (GS( K) ,K = 1 , NG ) 

WRITE*6,1008)(GSIK) ,K=1,NG) 

1008 FORMAT ( 6H G /10E12.5) 

WRITE( 6,996) 

READt 5 * 1 005 ) ( P AR P { K ) > K=1 * NG ) 

WR I TE ( 6,1009) t P ARP * K) , K=1 ,NG) 

ICO 9 FORMAT* 10E12.5) 

DO 18 K=1,NG 
M=M + 1 

V YY =EXP (R ZERO*P ARP ( K ) ) 

IF (YYY.GT.1.0E+27) GO TO 888 
IF *GS*K).GT.1.CE+10) GO TO 888 

18 G ( M ) = G 5 { K ) *Y YY 
GO TO 189 

688 G(M >=1.0E + 37 
189 CONTINtE 

WZERO=ALPHA*RU{ l)*RZERO**3 
DO 20 J=2 » N 

READt £ , 10C5) (PAR(K) ,K=1, NG ) 

WRITE* 6,1009) (PAR <K) ,K = 1,NG) 

DO 19 K = 1 , NG 
L =L + 1 

PtLfJ-1) =PARP ( K )~P AR * K) 

19 PARP*K)=PAR(K) 

Qt J-i )=ALPHA* IRO* J)-RQ( J-l ) ) 

20 L =L-NG 

DO 21 K = 1 » NG 
L=L+1 

21 P ( L t N ) =PAR ( K ) 

Q t N )=— £LP HA*RG ( N ) 

GO TO 100 
ENTRY CR 
D = 0 .0 

DO 32 K = 1,L 
BAD =0 * C 
DO 31 1 = 1 ,N 

31 BAD=BAOR( I)*P(K,I ) 

32 D=D+G( K )*EXP (-BAD) 

CON=ALCG( O/DBAR )+l • 0 
GO TO 100 

ENTRY fcR 
W=WZER0 
DO 41 1=1, N 
4 1 W=W+Ql I )*R(I )**3 
GO TO 100 
ENTRY EDIT 

WRITE* 6,1010) RZERC, ( R( I) ,1 = 1 ,N) 

1010 FORMAT *4HCRO=E12» 5/5H R { I )/< 1 OE 12. 5 ) ) 

D = 0 .0 
W=WZER C 
J =0 

M=M~N G 
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MP=M+I 
MPP =M / i +N G 
00 42 K*l,L 
BAD = 0 « C 
00 46 1 = 1 ,N 
46 BAD=BAC+R< 1I*P(K,I) 

A=G t K J^EXP (-BAD ) 

DO $( K l=A 

42 D = D+A 

uO 43 K=1»M,2 
J=J+1 

43 OOSU >=DOSUO+DOSIK+1) 

DO 44 K =M P , L 

J=J + 1 

44 OOSU ) =00 S ( K ) 

1011 FORMA T (7HCD0S(I)/(10E12*5i ) 

WRITE (6,1 Oil) (OOSU ) ,I=1,MPPJ 

DO 45 1 = 1, N 

45 W=W+Q( 1)«R(I)^3 
W=-K* l.OE + 5 

1012 FORMA T ( /3H W=E12,5,3H D=E12*5) 
WRITE(6,1012) W,D 

IOC RETURN 
END 


FORTRAN 63 COMPUTER CODE LISTING 


C0MM0N/MAX/N«Z<2)»KI 19)»Xt 2) *CDN«W 

DIMENSION RMI (20) ,RMA ( 20 ) , DER I 2U J *KOPT < 20 ) « RNGI 20 ) 

100 CALL INITIAL 

READ 101 ,EPS.nITMAX»ISW 
READ 1 02 » ( RM I ( I ) , I*1»N) 

PRINT102»(KMI ( I ) , I a l t N ) 

HE AO 102 * l KM A II > » I = 1 1 N I 
PRINT102 » IRMA ( I ),I=1,N) 

READ 102? (DERII ) , 1*1 . N ) 

PR INT102 * ( DEK I i ) * I “1 »N) 

101 FORMAT l fc 1 2 ••+ » 2 1 6 ) 

102 FORMAT (6612. 4» 

WMAX=-1 • E+200 $ NM1-N- 1 S DO 1 1*1 tMMl $R ( 1 I*RM1 1 I ) 

1 «NG( I )«RMHNI S. RID" R ( 1 ) -DERI 1 ) 

2 1*1 

3 R ( I ) =R ( I I+DERl I ) 

I F ( R ( I ) -RMA ( l ) 14,4,15 

4 IFIRIII-RI 1 + 1 1)5,5 1 15 

5 R ( N ) =RN1 =RNG I I ) $NI T=N I TMAX SCALL DR SC0N1=G0N-1 
IF (ABSFtCONl )-EPS)10,10,6 

6 R ( N ) = RN2 = RN'l +DER ( N ) iCALL DR iC0N2*C0N-l 
T r ("ABSF I C.nN2 ) “EPS ) 10 , 10 , 7 

7 R ( N ) =RN=RN 1-CON] *( RN1-RN2 ) / ( CON1-CON2 ) SCALL D.R SCON-CON-l 
IF (ABSF (CON ]“EPS)10,10,B 

8 NIT = N IT — 1 s IF(NIT)99,99»9 
99 STOP 

9 RN2 = RN1 5>C0N2=G0Nl $ RNi=RN $ CON1*CON $ GO 1 Q 7 

10 RN6 ( I ) =R ( N ) $ IF(R(N)-R(N-1))2»11,1L 

11 CALL WR $ IF ( I SW ) 12 » 1 2, 1 1 1 

111 PRINT 1I2,(W,(R(J),J=1,N)I 

112 FORMAT ( 10612*41 

12 I F ( W- WM A X ) 2 , 2 * 1 3 

13 W M A X = W $ DO 14 J*1,N 

14 ROPT ( J ) =R ( J ) $60 TO 2 

15 R ( I ) =RM I ( I ) $ lF(I-l)99,18*16 

16 I F ( R ( I )-R( 1-1 ) >17,18,18 

17 R( I )=R( 1-1 ) 

18 1=1 + 1 $ IFU-NJ 3,19,19 
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19 DU 20 1=1 *20 

20 R ( I ) = R0PT ( I > s CALL EDIT $G0 TO 100 

END 

SUBROUTINE INITIAL ~ 

COMM0N/MAX/N,MNOP»ALAM»R( 19 ) , DZZ Z ,CSO . CON , W 

DIMENSION R0( 19) ,GS(10) »Z ( 10) ,PAR( 10) ,G< 100) , PARPT10) * PC 100,19), 
10(19) f DOS ( 100 1 
' DIMENSION TITLE! 10) 

999 FORMAT ( 1H1 ) 

WRITE OUTPUT TAPE 6,999 

ALPHA= (4.0*3. 14159) /(453,5924E+5#3. ) 

1000 FORMAT ( 1 0A8 ) 

READ INPUT TAPE 5,1000, (TITLELJ ) ,1=1,10) 

WRITE OUTPUT TAPE 6 , 1000, ( T I TLfc ( I") , I * 1 , 10 ) 

1002 FORMAT I 3H0N= I 1 2 , 7H RZERO*E 1 2 . 5 ,6H 0BAR=E12.5) 

2002 FORMAT ( I12»2E12»5) 

READ INPUT TAPE 5 , 2002 , N, RZERO, DBAR 
WRITE OUTPUT TAPE 6 , 1 002, N * RZERO, D8AR 
L=M=KP=n 

DO 17 1=1 .N 

1003 FORMAT I 4HOKR=I5 »4H NG=15,4H R0«E12,5) 

READ INPUT TAPE 5,Z003*KR ,NG,ROKR 

2003 FORMAT 12I5»E12»5) 

KU(KK)=RUKR 

WRITE OUTPUT TAPE 6 , 1003, KR , NG,RO ( KR ) 

IF(KR-I) 901,1,901 
901 PRINT 901 

901 FORMAT! 31 H REGION NUMBERS NOT CONSECUTIVE) 

STOP 

1 READ INPUT TAPE 5 » 1 005 » < GS ! K ) » Z l K ) , K= 1 • NG) 

1004 FORMAT (5E 12.5 ) 

WRITE OUTPUT TAPE 6,998 
998 FORMAT! 4H GS) 

WRITE (6, 1004) ( GS (K ) *K=1 , NG ) 

WRITE OUTPUT TAPE 6,997 
997 FORMAT !4H Z) 

WRITE (6, 1004) ( Z ( K I ,K=1,NG) 

WRITE! 6,996 ) 

996 FORMAT ( 4H PAR) 

NG=2*NG 
DO 16 J=1,N 

READ INPUT TAPE 5 , 1 005 , ( P AR ( K ) , K= 1 , NG ) 

1005 FORMAT ( 1 0E8 .4 ) 

WRITF OUTPUT TAPE 6,1006, ! PAR < K > , K= l , NG ) 

1006 FORMAT! -10EI2.5) 

IF(J-l) 902,2,4 

902 PRINT 902 

902 FORMAT ( 1 6H J LESS THAN ONE) 

STOP 

2 00 3 K = ) ,NG » 2 
KP-KP+l 

3 GS!KP)«GS!KP)#EXPF(RZeRO*PAR(K) ) 

K P = 0 

4 IF(J-I) 5,7,12 

5 DO 6 K = 1 , NG , 2 

6 PAR ( K+l ) =PAR ( K ) *G0 TO 9 

7 DO 8 K = 1 , NG , 2 
M=M+ 1 $KP=KP+1 

A=G(M)=-GS!KP)fZ!KP)/!PAR!K )-PAR!K +1 I ) 
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M-M+l 

IF < 1—1 ) 902*903*904 

903 6 <M)=-A*EXPF (RZ£R 0 *IPAM K+ 1 )-PAR(K) > I SGO TO 8 

904 6!M)=-& 

8 CONTINUE 
KP=0 

9 IF(J-l) 902*10*14 

10 DO 11 K= 1 *N 6 

11 PARPl K ) “PAR ( K ) SGO TO 16 

12 DO 13 K = 1 *N 6 ♦ 2 

13 PAR ( K ) = PAk ( K+ ] ) 

14 DO 15 K= 1 ,NG 
L=L.+l 

PILvJ-1 ) sRAR.P f R ) “PAR C K ) 

15 PARPIKI-PARiK) $ L-l-NG 

16 CONTINUE 

DO 17 K* 1 *NG 
L=L +1 

17 P(t*N)=PAR(K) 

READ ( 5 * 2002 ) (NG) 

1007 FORMAT I AMONG- 1 10 } 

WR ITE ( 6 *1 007 ) (NG) 

READ! 5 * 1005 ) (GS IK ) *K= 1 »NG ) 

WRITE I 6 . 100 B HGS( K ) *K = i*NG) 

1008 FORMAT { 6 H 6 / 10 E 12 . 5 ) 

WRITE (6*996) 

READ! 5 * 1005 ) (PARP(K) *K=l*NGl 
WRITE 16 * 1 009 ) (PARP(K) *K= 1 *NG) 

[009 FORMAT I 10612 * 5 ) 

DO 18 K* 1 *NG 
M»M+i 

18 GIM)=GS(K)#EXPF»RZERQ*PARP(lU ) $WZERD=ALPHA*RO< 1 )*RZER 0**3 
DO 20 J = 2 *N 

READ! 5 ♦ 1005 )(PAR(R)*K*ltNG) 

WRITE ( 6 * 1009 ) ( PAR ( K ) *K= 1 *NG ) 

DO 19 K “ 1 » NG 3 . L*L *1 
P!L»J- 1 )=PARP(K)-PARIK) 

19 PAKP(K 1 *PAR(R) $ QU-l)s ALP HA* (ROT JJ-ROl J-l) ) 

20 l*L-NG 

DO 21 K-l *NG $ L=L +1 

21 P(L*N)*PAR(K ) S D(N)*-AtPMA*RQ<N) 

GO TO 100 

ENTRY OR 
D *0 *0 

DO 32 K-l * i. 

BAD®0®0 
DO 31 1 * 1 *N 

31 BAD*BAO+R(I)*P{K«n 

32 0*D+G(K )*EXPF(~8AD ) 

CON*LDGF( 0 /D 8 AR )+1 
GO TO 100 

ENTRY WR 

w*wzero 

DO 41 1 * 1 i N 

41 w*w+om*R( n *#3 
60 TO 100 
ENTRY EDIT 

WRITE OUTPUT TAPE 6* 1010*RZERG* (R (I ) t i*l » N) 
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1010 F0RMAT(4H0R0=E12*5/5H K( I l/(10E12.5> ) 

[)«0.0 S W=WZERO S J*0 S M=M-NG S MP»M+1 & MPP*M/2+NG 
DO 42 K = 1 * L 
BAD=0 o 0 
DO 46 1=1 »N 

46 RAD=BAD+R< I )#P<K, I ) 

fl=DOS ( K ) =G ( K ) #EXPF ( -BAD ) 

42 D=D+A 

DO 43 K=1 »M»2 
J = J + 1 

43 DOS ( J ) = DOS ( K 1 +DUS ( K+l ) 

DO 44 K*MP»L 

d = J + l 

44 DO S ( J ) = DO $ (K I 

1011 FORMAT ( 7H000S I I I/U0E12.5) ) 

WRITE OUTPUT < APE 6 ♦ 1 0 1 1 « < DOS ( I) » I * 1 » MP P } 

DO 45 1=1 .M 

45 W = W+t)(I)*R r )*#3 
W=-W*l .OE+5 

1012 F0RMAT(/3h W=E12.5 t 3H D=E12.5> 

WRITE OUTPUT TAPE 6,1012#W,D 

100 RETURN 
END 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, April 9, 1970, 

126-15. 
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